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Abstract. Optical flow cannot be computed locally, since only one independent measurement is available 

from the image sequence at a point, while the flow velocity has two components. A second constraint 
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is needed. A method for finding die optical flow pattern is presented which assumes that the apparent 
velocity of the brightness pattern varies smoothly almost everywhere in die image. An iterative implemen¬ 
tation is shown which successfully computes the optical flow for a number of synthetic image sequences. 
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included where the assumption of smoothness is violated at singular points or along lines in the image. 


Acknowledgement. This report describes research done at the Artificial Intelligence Laboratory of the 
Massachusetts Institute of Technology. Support for this research is provided in part by the Advanced 
Research Projects Agency of the Department of Defense under Office of Naval Research contract 
N00014-75-00643. 

© WfcSACHUSETTS INSTITUTE of technology 1980 



Determining Optical Flow 


Page / 


1. Introduction 

Optical flow is the distribution of apparent velocities of movement of brightness patterns in an image. 
Optical flow can arise from relative motion of objects and the viewer (Gibson 1950,1966). Consequently, 
optical flow can give important information about the spatial arrangement of the objects viewed and the 
rate of change of this arrangement (Gibson 1977). Discontinuities in the optical flow can help in segment¬ 
ing images into regions that correspond to different objects (Nakayama & Loomis 1974). Attempts have 
been made to perform such segmentation using differences between succcsive image frames (Jain et al. 
1977, Jain et al. 1979, Jain & Nagel 1979, Limb & Murphy 1975, Nagel 1977). Some recent papers have 
considered the problem of recovering the motions of objects relative to the viewer from the optical flow 
(Hadani et al. 1980, Kocndcrink & van Doom 1975 & 1976, Longuet-Higgins & Prazdny 1979, Prazdny 
1979 & 1980). In some cases information about the shape of an object may also be recovered (Kocndcrink 
& van Doom 1975 & 1976, Clocksin 1978). 

These papers begin by assuming that the optical flow has already been determined. Although some 
reference has been made to schemes for computing the flow from successive views of a scene (Fcnncma 
& Thompson 1979, Hadani et al. 1980), the specifics of a scheme for determining the flow from the image 
have not been described. Related work has been done in an attempt to formulate a model for the short 
range motion detection processes in human vision (Batali & Ullman 1980, Marr & Ullman 1979). 

The optical flow cannot be computed at a point in the image independently of neighboring points 
without introducing additional constraints, because the velocity field at each image point has two com¬ 
ponents while the change in image brightness at a point in the image plane due to motion yields only one 
constraint. Consider for example a patch of a pattern where brightness varies as a function of one image 
coordinate but not the other. Movement of the pattern in one direction alters the brightness at a particular 
point, but motion in the other direction yields no change. Thus components of movement in the latter 
direction cannot be computed locally. Additional constraints must be introduced to fully determine the 
flow. 


2. Relationship to Object Motion 

The relationship between the optical flow in the image plane and the velocities of objects in the three 
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dimensional world is not necessarily obvious. We perceive motion when a changing picture is projected 
onto a stationary screen, for example. Conversely, a moving object may give rise to a constant brightness 
pattern. Consider for example, a uniform sphere which exhibits shading because its surface elements are 
oriented in many different directions. Yet, when it is rotated, the optical flow is zero at all points in the 
image, since the shading does not move with the surface. Also, specular reflections move with a velocity 
characteristic of tine virtual image, not the surface in which light is reflected. 


3. The Restricted Problem Domain 

For convenience, we tackle a particularly simple world where the apparent velocity of brightness 
patterns can be directly identified with the movement of surfaces in the scene. To avoid variations in 
brightness due to shading effects we initially assume that the surface is flat. We further assume that 
the incident illumination is uniform across the surface. The brightness at a point in the image is then 
proportional to the reflectance of the surface at the corresponding point on the object. Also, we assume 
initially that reflectance varies smoothly and has no spatial discontinuities. This latter condition assures us 
that the image brightness is differentiable. 

In this simple situation, the motion of the brightness patterns in the image is directly determined by 
the motions of corresponding points on the surface of the object. Computing the velocities of points on 
the object is a matter of simple geometry once the optical flow is known. 


4. Constraints 


We will derive an equation that relates the change in image brightness at a point to the motion of the 
brightness pattern. Let the image brightness at the point ( x, y) in the image plane at time t be denoted by 
E(x, y, t). Now consider what happens when the pattern moves. The brightness of a particular point in 
the pattern is constant, so that 
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Using the chain rule for differentiation we see that, 

dE dx dE dy dE 
dx dt dy dt dt 

(See Appendix A for a more detailed derivation.) If we let 

dx 

u = — and v — 

d ft 


= 0 . 


dry 

dt’ 


then it is easy to see that we have a single linear equation in the two unknowns u and v. 


E x u -j- E y v -f ~ Et — 0, 

where we have also introduced the additional abbreviations E x , E y , and E t for the partial derivatives 
of image brightness with respect to x, y and t, respectively. The constraint on the local flow velocity 
expressed by this equation is illustrated in Figure 1. Writing the equation in still another way, 

(E x , E y ) • (u, v) = —E t . 

Thus the component of the optical flow in the direction of the brightness gradient (E x ,E y ) equals 

E t 

\J E l+ E i' 

We cannot, however determine the component of movement in the direction of the iso-brightness con¬ 
tours, at right angles to the brightness gradient. As a consequence, the flow velocity ( u, cannot be 
computed locally without introducing additional constraints. 


5. The Smoothness Constraint 

If every point of the brightness pattern can move independently, there is little hope of recovering the 
velocities. More commonly we view opaque objects of finite size undergoing rigid motion or deformation. 
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Figure 1. The basic rate of change of image brightness equation constrains the 
optical flow velocity. The velocity {u,v) has to lie along a certain line perpendicular 
to the brightness gradient vector ( E x ,E y ) in velocity space. 
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In this case neighboring points on the objects have similar velocities and the velocity field of the bright¬ 
ness patterns in the image varies smoothly almost everywhere. Discontinuities in flow can be expected 
where one object occludes another. 

One way to express the additional constraint is to limit the difference between the flow velocity at 
a point and the average velocity over a small neighborhood containing the point. Equivalently we can 
minimize the sum of the squares of the Laplacians of the x- and y-components of the flow. The Laplacians 
of u and v arc defined as 


V' 2 


_ d 2 u d 2 u 

U dx 2 dy 2 


and 


V 2 


_ d 2 v d 2 v 
dx 2 dy 2 


In simple situations, both Laplacians are zero. If the viewer translates parallel to a fiat object, rotates about 
a line perpendicular to the surface or travels orthogonally to the surface (assuming perspective projection), 
then die second partial derivatives of both u and v vanish. Note that our approach is in contrast with 
that of Fennema & Thompson (1979), who propose an algorithm that indirectly incorporates additional 
assumptions such as surface smoothness or object rigidity. 


6. Quantization and Noise 

Images may be sampled at intervals on a fixed grid of points. While tcssclations other than die 
obvious one have certain advantages (Mcrscrcau 1979, Gray 1971), for convenience we will assume that 
die image is sampled on a square grid at regular intervals. Let the measured brightness be Eij^ at the 
intersection of the i-th row and jf-tli column in the k -th image frame. Ideally, each measurement should 
be an average over the area of a picture cell and over the length of the time internal. In the experiments 
cited here we have taken samples at discrete points in space and time, instead. 

In addition to being quantized in space and time, the measurements will in practice be quantized in 
brightness as well. Further, noise will be apparent in measurements obtained in any real system. 
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7. Estimating the Partial Derivatives 

Wc must estimate the derivatives of brightness from the discrete set of image brightness measure¬ 
ments available. It is important that the estimates of E x , E y , and E t be consistent. That is, they should 
all refer to the same point in the image at the same time. While there arc many formulas for approximate 
differentiation (Conte & de Boor 1965, Hamming 1962) wc will use a set which gives us an estimate of E x , 
E y , E t at a point in the center of a cube formed by eight measurements. The relationship in space and 
time between these measurements is shown in Figure 2. Each of the estimates is the average of four first 
differences taken over adjacent measurements in the cube. 

E x Eij'k -f- Ei-\.Ei- 

-f-Ei.jH- i,M-i —Ei,j,k +1 —E i+l>ji k+i} 

Ey ^^{Ei+\,j,k — E itj)k + Ei+i t j+i t k — Eij+i'k 

+ ^*-F —Ei,j,k +1 +£*+i,i+i,M-1 Eij+-i t h-\-i} 

Et ^^{Ei,j,k-\~i Ei t j t k -\- Ei+i t j t k~\~i Ei-^\ j ik 

-{-■Eij+i'k+i — 

Here die unit of length is die grid spacing interval in each image frame and the unit of time is the image 
frame sampling period. 


8. Estimating the Laplacian of the Flow Velocities 

We also must approximate the Laplacians of u and v. One convenient approximation takes the 
following form 


V 2 u K{ui,j,ic — Uij'k) and V 2 v ^ k — 
where the local averages u and v are defined as follows 

Ui,j,k = g{ u i —lj,fc 4“ u i,j-\~\,k 4~ 4" u i,j— l,fc} 

4“ J24- 4- 4- «i+U—1,*} 

Vi,j,k =^{Vi-l,j,k 4- v i,j+l,k 4- «*+l ,j,k 4- v i,j-\,k} 

4- 4- Vi-i J+ i (fc 4- Wi4-i f i+i,A; 4- 
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Figu re 2. The three partial derivatives of image brightness at the center of the cube 
are each estimated from the average of first differences along four parallel edges 
of the cube. Here the column index j corresponds to the x direction in the image, 
the row index i to the y direction, while k lies in the time direction. 
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The proportionality factor k equals 3 if the average is computed as shown and assuming again that the 
unit of length equals the grid spacing interval. Figure 3 illustrates the assignment of weights to neighbor¬ 
ing points. 


9. Minimization 

The problem then is to minimize the sum of the errors in the equation for the rate of change of image 
brightness, 

Sf, = E x u -{- E y v -f- Et, 

and the estimate of the departure from smoothness in the velocity flow, 

e c = (« — u ) 2 + (t> — w) 2 . 

What should be the relative weight of these two factors? In practice the image brightness measurements 
will be corrupted by quantization error and noise so that we cannot expect & b to be identically zero. This 
quantity will tend to have an error magnitude that is proportional to the noise in the measurement. This 
fact guides us in choosing a suitable weighting factor, denoted by a 2 , as will be seen later. 

Let the total error to be minimized be 


g 2 = a 2g2 + g2. 

The minimization is to be accomplished by finding suitable values for the optical flow velocity (u, v). We 
differentiate 8 2 to obtain 

/5g2 

V- = - 2a 2 (u — u) + 2 (E x u + E y v + E t )E x 
ou 

/9g2 

^ = ~ 2a\v - v) + 2 (E x u + E y v + E t )E y . 
ov 

Setting these two derivatives equal to zero leads to two equations in u and v. 
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(a 2 + El)u + E t E y v = (a 2 u - E x E t ) 
E z E v u + (a 2 + E 2 )v = (a 2 6 - EyE t ). 
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Figu re 3. The Laplacian is estimated by subtracting the value at a point from a 
weighted average of the values at neighboring points. Shown here are suitable 
weights by which values can be multiplied. 
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The determinant of the coefficient matrix equals q 2 (q 2 -f E\ + £ 2 ). Solving for u and v we find that 

(q 2 -j- E 2 -f- E 2 )u == -j- (<a 2 -{- E 2 )u ExE y v E x Et 
(a 2 + £ 2 + £> = - + (a 2 + jE*)i> - EyE t . 


10. Difference of Flow at a Point from Local Average 

These equations can be written in the alternate form 

(a 2 + E\ + E\\u - u) = - E x [E x u + E y v + E t \ 

(a 2 + E 2 X + El)(v - v) = - + 3]• 

This shows that the value of the flow velocity ( u, v) which minimizes die error 6 2 lies in the direction 
towards the constraint line along a line that intersects the constraint line at right angles. This relationship 
is illustrated geometrically in Figure 4. The distance from the local average is proportional to the error 
in the basic formula for rate of change of brightness when ii, v arc substituted for u and v. Finally we 
can see that a 2 plays a significant role only for areas where the brightness gradient is small, preventing 
haphazard adjustments to die estimated flow velocity occasioned by noise in the estimated derivatives. 
This parameter should be roughly equal to die expected noise in the estimate ofi? 2 -f- E 2 . 


I t. Constrained Minimization 

When we allow a 2 to tend to zero we obtain the solution to a constrained minimization problem. 
Applying the method of Lagrange multipliers (Russell 1976, Yourgau & Mandelstam 1968) to the 
problem of minimizing 8 2 while maintaining gj, = 0 leads to 

(E\ -f E'l)(u - u) = - E x \E x u + E y v + E t ] 

{E 2 X + E 2 y ){v — v) = — E y [E x ii + E y v + E t \. 

We will not use tlicsc equations since we do expect errors in the estimation of the gradient (E X) E y ). 
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12. Iterative Solution 

We now have a pair of equations for each point in the image. It would be very costly to solve these 
equations simultaneously by one of the standard methods, such as Gauss-Jordan elimination (Hamming 
1962, Hildebrand 1952). The corresponding matrix is sparse and very large since the number of rows and 
columns equals twice the number of picture cells in the image. Iterative methods, such as the Gauss-Seidel 
method (Hamming 1962, Hildebrand 1956), suggest themselves. We can compute a new set of velocity 
estimates {u n ^ \ i> n-t_1 ) from the estimated derivatives and the average of the previous velocity estimates 
(u n , v n ) by 

„»+< —u n - E x {E x u n + E y v" + E t ]/(a 2 + B 2 + E 2 y ) 

„»+i =0- - E y [E x u n + E y v n + E t ]/(a 2 + E\ + E 2 ). 

It is interesting that the new estimates at a particular point do not depend directly on the previous 
estimates at the same point. 


13. Filling In Uniform Regions 

In parts of the image where the brightness gradient is zero, the velocity estimates will simply be 
averages of the neighboring velocity estimates. There is no local information to constrain the apparent 
velocity of motion of the brightness pattern in these areas. Eventually the values around such a region will 
propagate inwards. If the velocities on the border of the region arc all the same, then points in the region 
will be assigned that value too after a sufficient number of iterations. Velocity information is thus filled in 
from the boundary of a region of constant brightness. 


14. Number of Iterations 

If the values on the border are not all the same, the values filled in will correspond to the solution 
of the Laplace equation for the given boundary condition (Ames 1977, Milne 1953, Richtmycr& Morton 
1957). The progress of this filling in phenomena is similar to the propagation effects in the solution of the 
heat equation, where the time rate of change of temperature is proportional to the I aplacian. This gives 
us a means of understanding the iterative method in physical terms and of estimating the number of steps 
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required. The number of iterations should be larger than the cross-section of the biggest region that must 
be filled in. If the size of such regions is not known in advance one may use the cross-section of the whole 
image as a conservative estimate. 


15. Tightness of Constraint 

When brightness in a region is a linear function of the image coordinates we can only obtain the 
component of optical flow in the direction of the gradient. The component at right angles is filled in from 
the boundary of the region as described before. In general the solution is most accurately determined 
in regions where die brightness gradient is not too small and varies in direction from point to point. 
Information which constrains both components of the optical flow velocity is then available in a relatively 
small neighborhood of each point. I'oo violent fluctuations in brightness on the other hand are not 
desirable since the estimates of derivatives will be corrupted due to undersampling and aliasing. 


16. Choice of Iterative Scheme 

As a practical matter one has a choice of how the iterations are to be interlaced with the time steps. 
On the one hand, one could iterate until the solution has stabilized before advancing to the next image 
frame. On the other hand, given a good initial guess one may need only one iteration per time-step. A 
good initial guess for the optical flow velocities is usually available from the previous time-step. 

The advantages of the latter approach include an ability to deal with more images per unit time and 
better estimates of optical flow velocities in certain regions. Areas in which the brightness gradient is small 
lead to uncertain, noisy estimates obtained partly by filling in from the surroundings. These estimates are 
improved by considering further images. The noise in measurements of the images will be independent 
and tend to cancel out. Perhaps more importantly, different parts of the pattern will drift by a given point 
in the image 

A practical implementation would most likely employ one iteration per time step for these reasons. 
We illustrate both approaches in the experiments. 
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17. Experiments 

The iterative scheme has been implemented and applied to image sequences corresponding to a 
number of simple flow patterns. The results shown here are on a relatively small image of 32 by 32 
picture cells. The brightness measurements were intentionally corrupted by approximately 1% noise and 
then quantized into 256 levels to simulate a real imaging situation. The underlying surface reflectance 
pattern was a linear combination of spatially orthogonal sinusoids. Their wavelength was chosen to give 
reasonably strong brightness gradients without leading to undersampling problems. Discontinuities were 
avoided to ensure that the required derivatives exist everywhere. 

Shown in Figure 5 for example arc four frames of a sequence of images depicting a sphere rotating 
about an axis inclined towards the viewer. A smoothly varying reflectance pattern is painted on the surface 
of the sphere. The sphere is illuminated uniformly from all directions so that there is no shading. 


18. Results 

The first flow to be investigated was a simple linear translation of the entire brightness pattern. Hie 
resulting computed flow is shown as a needle diagram in Figure 6 for 1, 4, 16, and 64 iterations. The 
estimated flow velocities are depicted as short lines, showing the apparent displacement during one time 
step. In this example a single time step was taken so that the computations are based on just two images. 
Initially the estimates of flow velocity are zero. Consequently the first iteration shows vectors in the 
direction of the brightness gradient. Later, the estimates approach the correct values in all parts of the 
image. Few changes occur after 32 iterations when the velocity vectors have errors of about. 10%. The 
estimates tend to be too small, rather than too large, perhaps because of a tendency to underestimate the 
derivatives. The worst errors occur, as one might expect, where the brightness gradient is small. 

In the second experiment one iteration was used per time step on the same linear translation image 
sequence. The resulting computed How is shown in Figure 7 for 1, 4, 16, and 64 time steps. The estimates 
approach the correct values more rapidly and do not have a tendency to be too small, as in the previous 
experiment. Few changes occur after 16 iterations when the velocity vectors have errors of about 7%. The 
worst errors occur, as one might expect, where the noise in recent measurements of brightness was worst. 
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A 


C D 

Figure 5. Four frames out of a sequence of images of a sphere rotating about an 
axis inclined towards the viewer. The sphere is covered with a pattern which 
varies smoothly from place to place. The sphere is portrayed against a fixed, lightly 
textured background. Image sequences like these are processed by the optical flow 
algorithm. 
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Figure 6. Flow pattern computed for simple translation of a brightness pattern. The 
estimates after 1, 4, 16, and 64 iterations are shown. The velocity is 0.5 picture 
cells in the x direction and 1.0 picture cells in the y direction per time interval. 
Two images are used as input, depicting the situation at two times separated by 
one time interval. 


A. /. Memo 572 


Horn & Schunck 
















Determining Optical Flow 


Page 17 


) >. \ 

i'l'k'kN.'.'k'k’l, 

UuuU 

\\ww\W 

iif 

77/////// 


J > \ \ w \ \ \ 
\\wwwV 

• j ■>.*«. ^ k k & • 

y//M 


^•«ftftftfta« 

\\'k>.'k'k'k\\ 

IWwwW 

\'l'k\'k'k'k\\ 



*, *,4 4444 4 i, 

-. Ilss///l1+ 

* s r-.-'r " w'w' 

l\\\k'»k'.'k'i 

Uuu\ 

Uu\u\ 

l^'k'kSk'.\ V k \ 
*,'» ■»»» \ 

mm 

mr/rni 

11/mm//ll 

/ —-.—— 

\ > -.. \ 

[Uu'.wU 

\uvunl 

UunswV\ 

y .. 

Uwwt 
WM 


'»■•'» • ..'k'k 

\WkkN.\\ 
V'k'kk^N.W, 

l'k\\.\'k'k\. 
'k'k v »k-»'*k\l 

'^V/V/Aii 

mm 



i 


mm 
yjmfa 

.> 

U ^ ^ ; u 

U H \ u \ I 

^us\Ul 



j //////i i 

mm, 

it It It ill 
i m i ill/. 

\ \ V k k k l i ill 

WM ix 
if 

‘■mm 

1 1 ill U 11 11 
i 1 1 M M 1 1 1 1 
UvuillJ/ 

- V V v k I » l 1 ill 

\\ j p//// / /j 1 l l pl 



fl| 

I U n m 

\\UUll. 

\, k i * » t i i i 

ifiPPt 

jpjM. 


i \\\\ww\ 
Uukun 
* * * i i. 


| u # 

4 *, *m 4 <r 4 4,4, 

IP 

PJM 

n v; 


hiii 


\'VZ/////i 

i 1 i i i i i 1 


1 l l v l 1 l 
l, k * * * * J J I, 

U/SSS/J.% 
!/////1,1,1 

m/rn/z/A 


*.*.* 44444 4 

m 

Uu u i i 


iMmmM 


iHiPity/iPP/YPAPPPP/PPPiPJA 

PAAAPAAPAAPAAAAAPAPAAPAA 

FitFiFl nil ii nil 1 1 11111 in y / j i 

llJlililillililTlTlilli 111 11 111 1 

wm PAA PA PA PPAPAtfAPA A 
mmPPAtPmPAPAPAPAi 

J iiiifii ill JiiiiAiiiiiiiiiiiAiii 

AiAAiAiAiAiAiAiAAAiAiAAAAiAiAAi 

m PPAiPAPAiPAPAPAiiPA 

////////////////////////7777777 

WMttWMWMWMWAPA 

Jll.J!liJ/illlllllTll/llJJJl//f/ 

J.l 111J11 /// itititijjiiiiiiiiiii 
11 ill ill ill iiiiivjitjitiiiiit/1 

riJirriririririririViririririri 


YiVPPYPPiWPiPPYPiiPPj 

PAPAPAPAPAi APAPAAAPA 

vjiWtWiWiWiWiWiYiWDYj 

HttfflAAAPAPAAAPAPPAA 

PIWIWIPPPPPPPPPPPPAPPPPPPPPA 

AAAAiAAAAiWiAii AiAAiAiAiA 

wmPPA PPA PA PAM A PAP 

Jill Ji il ill 1 j ii Till iiiiiiiii 111 
11111 Till 1111111111llllllllllll 
ill 111111111111 liiiiiilTii riiii 

11111111111111111 ii ii Ai ii ii AiAi 
Til 1111111111lliiil Til TfWrrn V 
1111 il 111 ill/ill il iiiiiiiiiiiii 
111111111 111 111 1111 Till 7777777V 
/ /// 111J11111lllTJlTlll 11111111 
ill 11111 1111 il 1 li ii z/iiiiiiiiii 

ttfflttMAAPAtAAPAPAA 

mmPAPtPAPPPAPAPAPAP, 


Figure 7. Flow pattern computed for simple translation of a brightness pattern. The 
estimates after 1, 4, 16, and 64 time steps are shown. Here one iteration is used 
per time step. Convergence is more rapid and the velocities are estimated more 
accurately. 
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While individual estimates of velocity may not be very accurate, the average over the whole image was 
within 1% of the correct value. 

Next, the method was applied to simple rotation and simple contraction of the brightness pattern. 
The results after 32 time steps are shown in Figure 8. Note that the magnitude of the velocity is propor¬ 
tional to the distance from the origin of the flow in both these cases. 

In all of the examples so far the Laplacian of both flow velocity components is zero everywhere. We 
also studied more difficult cases where this was not the case. In particular, if we let the magnitude of the 
velocity vary as the inverse of the distance from the origin we generate flow around a line vortex and two 
dimensional flow into a sink. The computed flow patterns arc shown in Figure 9. Here the computation 
involved many iterations based on a single time step. The worst errors occur near die singularity at the 
origin of flic flow pattern, where velocities are found which arc much larger than one picture cell per time 
step. 

Finally we considered rigid body motions. Shown in Figure 10 are the flows computed for a cylinder 
rotating about its axis and for a rotating sphere. In both cases the Laplacian of the flow is not zero and in 
fact the Laplacian of one of the velocity components becomes infinite on flic occluding bound. Since the 
velocities themselves remain finite, reasonable solutions are still obtained. The correct flow patterns are 
shown in Figure 11. Comparing the computed and exact values shows that the worst errors occur on the 
occluding boundary. These boundaries constitute a one dimensional subset of the plane and so one can 
expect that the relative number of points at which the estimated flow is seriously in error will decrease as 
the resolution of the image is made finer. 

In Appendix B it is shown that there is direct relationship between the Laplacian of flic flow velocity 
components and flic Laplacian of the surface height. This can be used to see how our smoothness con¬ 
straint will fair for different objects. For example, a rotating polyhedron will give rise to a flow which has 
zero Laplacian except on the image lines which arc the projections of flic edges of the body. 


19. Summary 

A method was developed for computing optical flow from a sequence of images. It is based on the 
observation that the flow velocity has two components and that flic basic equation for the rate of change 
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Figure 8. Flow patterns computed for simple rotation and simple contraction of a 
brightness pattern. In the first case, the pattern is rotated about 2.8 degrees per 
time step, while it is contracted about 5% per time step in the second case. The 
estimates after 32 time steps are shown. 
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Figure 9. Flow patterns computed for flow around a line vortex and two dimensional 
flow into a sink. In each case the estimates after 32 iterations are shown. 
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Figure 10. Flow patterns computed for a cylinder rotating about its axis and for a 
rotating sphere. The axis of the cylinder is inclined 30 degrees towards the viewer 
and that of the sphere 45 degrees. Both are rotating at about 5 degrees per time 
step. The estimates shown are obtained after 32 time steps. 
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Figure 11. Exact flow patterns for the cylinder and the sphere. 
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of image brightness provides only one constraint. Smoothness of the flow was introduced as a second 
constraint. An iterative method for solving the resulting equation was then developed. A simple im¬ 
plementation provided visual confirmation of convergence of the solution in the form of needle diagrams. 
Examples of several different types of optical flow patterns were studied. These included cases where the 
Laplacian of the flow was zero and cases where it became infinite at singular points or along bounding 
curves. 
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Appendix A. Rate of Change of Image Brightness 


Consider a patch of the brightness pattern that is displaced a distance 8x in the ^-direction and 8y in 
tlte y-direction in time St. The brightness of the patch is assumed to remain constant so that 

E(x, ij, t) = E(x -{-Sx,y-\- 8y, t -\- St). 


Expanding the right hand side about the point (x, y, t) we get, 

. lS t . dE dE . s dE . 

E{x, y, t ) = E(x, y, t) -f 8x— + 6y-~ + St — + e. 

Where e contains second and higher order terms in 8x, 8y, and 8t. After subtracting E(x, y, t) from both 
sides and dividing through by 8t we have 


SxdE 8ydE dE 
6t dx 8t dy dt 


+ 0 (St) = 0, 


where 0(<5£) is a term of order 8t, and we assume that 8x and 8y vary as 6t. In the limit as 8t —> 0 this 
becomes 

dE dx dE dy dE ^ 

dx dt dy dt dt 
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Appendix B. Smoothness of Flow for Rigid Body Motions 


Let a rigid body rotate about an axis ( lo x , lo v , lo z ), where the magnitude of die vector equals the 
angular velocity of the motion. If this axis passes through the origin, then the velocity of a point ( x , y, z) 
equals the cross product of (to x , o)y, oj z ), and ( x , y, z). There is a direct relationship between the image 
coordinates and the x and y coordinates here if we assume diat the image is generated by orthographic 
projection. The x and y components of the velocity can be written, 

u —u)yZ — l o z y 
v =U) Z X — uo x z. 

Consequently, 

V 2 u = -{- u) y V 2 z 
V 2 n = — u) x V 2 z. 

I his illustrates diat the smoothness of the optical flow is related directly to die smoothness of the rotaung 
body and that die Laplacian of the flow velocity will become infinite on the occluding bound, since the 
partial derivatives of z with respect to x and y become infinite diere. 
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